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A scheme is presented and software is 
documented for representing as integers 
input decimal numbers that have been 
stored in a computer as double precision 
floating point numbers and for carrying 
out multiplications, additions and subtrac- 
tions based on these numbers in an exact 
manner. The input decimal numbers must 
not have more than nine digits to the left 
of the decimal point. The decimal fractions 
of their floating point representations are 
all first rounded off at a prespecified loca- 
tion, a location no more than nine digits 
away from the decimal point. The number 
of digits to the left of the decimal point for 
each input number besides not being 
allowed to exceed nine must then be such 



that the total number of digits from the 
leftmost digit of the number to the location 
where round-off is to occur does not 
exceed fourteen. 
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1. Introduction 

Mainstream computers base integer and floating 
point arithmetic on fixed word lengths. As a conse- 
quence, only values with a limited number of signifi- 
cant digits can be represented directly so that the 
results of arithmetic operations may have to be round- 
ed off or truncated. Such errors can be avoided or, at 
least, mitigated, by implementing special algorithms 
for the execution of arithmetic operations, A fully 
"exact arithmetic", however, would have to be based on 
quotients of integers for representing numerical values. 
In any case, a final limitation is due to finite memory. 

The need for exact arithmetic became apparent dur- 
ing the development of software for generating triangu- 
lar and tetrahedral nets from very large point sets. 
Typically, this need is not due to high accuracy require- 
ments for results — the input data are often noisy or 



given up to only a few significant digits — ^but is rather 
due to the need to maintain the consistency of a combi- 
natorial structure. Building or manipulating such geom- 
etry-based combinatorial structures requires the calcu- 
lation of indicators such as determinants in order to 
evaluate their sign and to check for zero values: round- 
off may lead to a false sign or zero value. An example 
considered later in this paper is to decide whether four 
given spatial points are coplanar. The approach of using 
exact computations for implementing computational 
geometry algorithms in a robust manner has been 
addressed in [2], [3], [4], [5], [6], [7]. Some computa- 
tional geometry implementations ([2], [3], [4], [5]) 
reduce computational effort by utilizing exact arith- 
metic selectively whenever a decision might be affect- 
ed by round-off. 

In this paper, we document software for exact integer 
arithmetic, accommodating an indeterminate number of 
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digits, for multiplication, addition, subtraction, but 
excluding division. We found that, in many computa- 
tional geometry applications, decision variables such as 
determinants can be calculated without division. Also 
the sign of a decision variable stated as a quotient but 
not evaluated is readily derived from the signs of the 
numerator and denominator. 

We also describe a preprocessing step, called "two- 
integer decomposition", which leads from floating 
point input to one composed of integers only. At the 
root of this step lies the concept of space as an integer 
grid of points, all of which have integer coordinates in 
some shared unit. After completing the transition from 
floating point numbers to intermediate representations 
as pairs of integers — ^prompted by the fact that Fortran 
77 does not provide a double precision integer 
format — a polynomial decomposition creates the num- 
ber representations to be used in the exact arithmetic 
calculations. Software for this preprocessing step 
together with software for exact integer arithmetic has 
been successfiilly incorporated into several computa- 
tional geometry related programs such as REGTET [1]. 

In what follows, a "standard computer" is a comput- 
er that uses 64 bits of storage for a double precision 
number and 32 bits for an integer. Given a standard 
computer, even though it may not store exactly an input 
decimal number as a double precision floating point 
number, it is safe to assume that the number will be rep- 
resented as accurately as possible by a double precision 
floating point number up to its fourteenth significant 
digit. 



2. Two-Integer Decomposition 

Let x(/), / = 1, ..., n, be a double precision array into 
which input numbers Xy, / = 1, ..., n, have been read. The 
two-integer decomposition process is a preprocessing 
step that takes place before any computations based on 
the input data are carried out. It rounds off the numbers 
in the array at a prespecified location of their decimal 
fractions and decomposes each rounded off number 
into two integers that are saved in integer arrays, say 
/x(/), ix2{i), i = 1, ..., 77. The rounded off numbers are 
then saved in array x(/), / = 1, ..., n. 

Given integers k, I, I < k < 9, < I < 9, k + I < 14, 
and assuming each input number x„ / = 1, ..., n, has no 
more that k digits to the left of the decimal point, each 
number x(/), / = 1, ..., n, is rounded off at the /th digit of 
its decimal fraction and decomposed into two integers 
in one of two ways according to its size. If the 
absolute value of x(/) times (lO.Oc/Oy is less than 



2^^1073741824), x(0 is multiphed by (10.0^0/ and 
rounded off at the decimal point. The resulting integer 
is then placed in /x(/) while /x2(/) is set to zero. Finally, 
x(/) is redefined to be the double precision value of 
integer /x(/) divided by (10.0<iO)^. On the other hand, if 
the absolute value of x(/) times (10.0<iOy exceeds or 
equals 2^^, x(i) is truncated at the decimal point. The 
resulting integer (absolute value less than 2^^ since k < 
9) is placed in /x(/). In addition, the signed decimal 
fraction obtained by subtracting the double precision 
value of this integer from the initial value of x(/) is mul- 
tiplied by (10.0<iO)^ and rounded off at the decimal 
point. The resulting integer (absolute value less than 2^^ 
since / < 9) is placed in ix2(i). Next, x(/) is redefined to 
be the double precision value of integer ;x(/) plus the 
value obtained by dividing the double precision value 
of integer /x2(/) by (lO.Oc/Oy. Finally, if the integer 
ix2(i) is zero then /x2(/) is set to 2^^ so that ix2(i) is zero 
if and only if the initial absolute value of x(/) (before 
the two-integer decomposition process) times (10.0<iO)^ 
is less than 2^^. 

The following is Fortran code for carrying out the 
two-integer decomposition process. Variables are either 
integer or double precision following convention. It is 
noted that while some loss in precision may occur at the 
time the input numbers are read and transformed into 
double precision floating point numbers, some addi- 
tional loss in precision may occur here as well when the 
decimal point in a number is shifted by dividing or mul- 
tiplying it by a multiple of 10.0<iO, when the signed 
decimal fraction of a number is obtained by truncating 
the number at its decimal point and subtracting the 
result from the initial value of the number, and when a 
number is rounded off with the two-integer decomposi- 
tion process. However once the two-integer decompo- 
sition process is completed all computations that fol- 
low, exact and otherwise, are carried out in terms of the 
arrays x(/), /x(/), /x2(/), / = 1, ..., n, under the assump- 
tion that for the purposes of the user for each 
/,/=!,...,«, x(/) represents closely enough the input 
number X/ rounded off at the /th digit of its decimal frac- 
tion, and an integer (not necessarily stored by the com- 
puter) in terms of /x(/) and /x2(/) represents closely 
enough the input number x^ times 1 0^ rounded off at the 
decimal point. 

mfull=1073741824 
if(Ut.O .or. l.gt.9) stop 10 
isclxi = 1 
dscle= l.OdO 
if(l.eq.O) go to 200 
do 100 i= 1,1 
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isclu= 10*isclu 
dscle= 10.0dO*dscle 
100 continue 
200 continue 

dfull = dble(mfull) 
dfill=dfull/dscle 
do300i= l,n 
ix2(i) = 

if(dabs(x(i)).lt.dfill) then 
ix(i) = idnint(dscle*x(i)) 
if(iabs(ix(i)).lt.mfull) then 
x(i) = dble(ix(i))/dscle 
go to 300 
endif 
endif 

if(dabs(x(i)).ge.dfull) stop 20 
ix(i) = idint(x(i)) 
if(iabs(ix(i)).ge.mfull) stop 30 
decml = (x(i) - dint(x(i)))*dscle 
ix2(i) = idnint(decml) 
if(iabs(ix2(i)).eq.O) then 
x(i) = dble(ix(i)) 
ix2(i) = mfull 
else 

x(i) = dble(ix(i)) + (dble(ix2(i))/dscle) 
endif 
300 continue 



3. Polynomial Decomposition 



Given an integer /, < / < 9, let Xy, / = 1, ..., n, be 
input numbers whose double precision floating point 
representations have been rounded off at the /th digit of 
their decimal fractions through the two-integer decom- 
position process. Let x{i\ ix(i), ix2(i), / = 1, ..., n be the 
arrays produced by the two-integer decomposition 
process that contain the rounded off numbers and the 
two-integer decompositions. For each /, i = 1, ..., n, an 
integer J(/, /) is symbolically defined as follows (its 
actual value is not necessarily computed or stored by 
the computer). If ix2(i) equals zero then J(i, I) is set 
equal to ix{i). If ix2{i) equals 2^^ then J(z, I) is set equal 
to ix{i) • 10^. Finally if a2(/) is neither zero nor 2^^ then 
J(i, I) is set to ix{i) • 10^ + ix2(i). In all cases for each /, 
/ = 1, ..., 77, J(i, I) is considered to approximate closely 
enough (for the purposes of the user) the input number 
Xy times 10^ rounded off at the decimal point. 

Set M to 2^^ Given J(z, 1), \ <i <n, the polynomial 
decomposition process is a procedure (presented below 
in the form of Fortran subroutine decmp2) that decom- 
poses the integer J(/, /) into a unique collection of inte- 



gers isga, isga in {-1,0, 1 } , ika, ika > 0, a^, < a^ < M, 
k= \, ..., ika, such that J(/, I) equals isga^L'^i a^ -A/^"^), 
isga the sign of J(z, I), Integers a^, k = 1, ..., ika, are 
saved in an integer array, say ia(k), k = I, ..., ika, and 
the collection of integers isga, ika, ia(k), k= I, ..., ika, 
and the symbolic expression isga(l!^i ia(k) • 71/"^) are 
then called, respectively the polynomial decomposi- 
tion and the symbolic polynomial representation of 
J(i, I), with isga as the sign of the representation. For 
each /, / = 1, ..., n, the polynomial decomposition of the 
integer J(i, I) is identified each time an exact computa- 
tion involving additions, subtractions, or multiplica- 
tions is required that references the input number x,. 
During one such computation, for each /, 1 < i < n, if 
the number x, is referenced in the computation, once the 
polynomial decomposition of the corresponding integer 
J(i, I) is identified, each reference of x, in the computa- 
tion is replaced by the symbolic polynomial representa- 
tion of J(i, I). The computation then takes effect as a 
sequence of additions, subtractions, or multiplications 
of symbolic polynomial representations with the final 
result being itself the symbolic polynomial representa- 
tion of some integer. This final result can usually be 
used in only one of two ways. If it is known that for 
some positive integer m the integer that is equal to the 
final symbolic polynomial representation is approxi- 
mately equal to the product of (10^)'" and the true value 
of the computation, then this integer is computed 
approximately as a double precision floating point 
number from its symbolic polynomial representation 
and the true value of the computation is then approxi- 
mately obtained by dividing it by {(lO.OdO))"'. On the 
other hand, if the purpose of the computation is simply 
that of obtaining the sign of the true result then the sign 
of the final symbolic polynomial representation is a sat- 
isfactory answer. 

The concepts of polynomial decomposition and sym- 
bolic polynomial representation defined above for 
J{i, I), 1 <i<n, can also be defined for any integer K 
(not necessarily stored by the computer) in the same 
manner. Accordingly, the following is a Fortran subrou- 
tine called decomp for finding the polynomial decom- 
position isga, ia{k), k= 1,2, (ika is already known to 
equal 2) of an integer iwi (stored by the computer) with 
absolute value less than 2^^. Here mhalf equals 
2^'(=32768). 

subroutine decomp(ia, isga, iwi, mhalf) 
integer ia(*), isga, iwi, mhalf, ivi 
if(iwi.gt.O) then 

isga = 1 

ivi = iwi 
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elseif(iwi.lt.O) then 

isga=-l 

ivi = -iwi 
else 

isga = 

ia(l) = 

ia(2) = 

return 
endif 

ia(2) = ivi/mhalf 
ia(l) = ivi - ia(2)*mhalf 
return 
end 

In particular if isclu is set to 10^ then isclu is less than 
2^^ (since / < 9) so that the polynomial decomposition 
isgu (equal to 1 ), iu(i), i= 1,2, (iku is already known to 
equal 2) of isclu can be obtained by calling subroutine 
decomp with a Fortran instruction as follows. 

call decomp(iu, isgu, isclu, mhalf) 

Finally, the following is a Fortran subroutine called 
decmp2 for finding the polynomial decomposition isga, 
ika, ia(k), k= I, .,., ika, of the integer J(z, I), \ <i <n. 
Here iwi equals ix(i), iwil equals ix2(i), mhalf equals 
1}^, mfull equals 2^^, and /w(A:), ^ = 1, 2, is an array such 
that the polynomial decomposition of 10^ is /w(l), m(2) 
(isgu and iku are already known to equal 1 and 2, 
respectively). In addition, it is assumed that subroutines 
mulmul and muldif (presented below) exist for multi- 
plying and subtracting, respectively, two symbolic 
polynomial representations, 

subroutine decmp2(ia, isga, ika, iwi, iwi2, mhalf, mfull, iu) 

integer nkmax 

parameter (nkmax=5) 

integer ia(*), isga, ika, iwi, iwi2, mhalf, mfull, iu(*) 

integer ie(nkmax), io(nkmax), isge, isgo, ike, iko, isgu, iku 

call decomp(ia, isga, iwi, mhalf) 

ika = 2 

if(iwi2.ne.O) then 

isgu= 1 

iku = 2 

call mulmul(ia, iu, ie, isga, isgu, isge, ika, iku, 

* ike,nkmax,mhalf) 
if(iwi2.eq. mfull) iwi2 = 

call decomp(io, isgo, iwi2, mhalf) 

isgo = -isgo 

iko = 2 

call muldif(ie, io, ia, isge, isgo, isga, ike, iko, ika, 

* nkmax,nihalf) 



endif 
return 
end 



4, Multiplying Symbolic Polynomial 
Representations 

Given the polynomial decompositions isga, ika, 
ia(k), k= I, ,.., ika, and isgb, ikb, ib(k), k= 1, ,.,, ikb, of 
two integers K^ and K2, respectively, the following is a 
Fortran subroutine called mulmul that produces the 
polynomial decomposition isgo, iko, io(k), k = 1, .,., 
iko, of the integer K^ • K2 by multiplying the symbolic 
polynomial representation ofK^ by that of A^2 (as poly- 
nomials) to produce a symbolic polynomial representa- 
tion of K^ • K2 from which the polynomial decomposi- 
tion of Ky • K2 can be obtained. Here nkmax is the 
dimension of the arrays ia, ib, io in the calling routine 
and mhalf cqu^ih l}^. It is noted that the value of mhalf 
is of importance here since given integers i, j, \ <i < 
ika, \<j< ikb, then < ia{i) < 2^\ < ibij) < 2^\ so 
that the product iaif) • ib(j) is less than 2^^ and therefore 
can be stored in a 32 bit integer word, 

subroutine mulmul(ia, ib, io, isga, isgb, isgo, ika, ikb, iko, 
* nkmax,mhalf) 

integer ia(*), ib(*), io(*) 

integer isga, isgb, isgo, ika, ikb, iko, nkmax, mhalf 
integer i, ipt, ipr, ikol, k, j 
if(isga.eq.0.or.isgb.eq.O)then 

isgo=0 

iko = 2 

io(l) = 

io(2) = 

return 
endif 

iko = ika + ikb 
if(iko.gt. nkmax) stop 110 
if(isga.gt.O)then 

if(isgb.gt.O)then 
isgo = 1 

else 

isgo =-1 

endif 
else 

if(isgb.gt.O)then 
isgo =-1 

else 

isgo = 1 

endif 
endif 
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ikol =iko - 1 




ipr = 




do 2001=1, ikol 




ipt = ipr 




k = i 




do 180 j = l,ikb 




if(k.lt. l)goto 190 




if(k .gt. ika) go to 150 




ipt = ipt + ia(k)*ibG) 


150 


continue 




k = k- 1 


180 


continue 


190 


continue 




ipr = ipt/mhalf 




io(i) = ipt - ipr*mhalf 


200 


continue 




io(iko) = ipr 




if(ipr.ge.mhalf) stop 120 




ikol =iko 




do300i = ikol,ika+l,-l 




if(io(i) .ne. 0) go to 400 




iko = iko - 1 


300 


continue 


400 


continue 




return 




end 



5, Subtracting Symbolic Polynomial 
Representations 

Given the polynomial decompositions isga, ika, 
ia{k), k=\, ..., ika, and isgb, ikb, ib{k), k=\, ..., ikb, of 
two integers K^ and K2, respectively, the following is a 
Fortran subroutine called muldif that produces the 
polynomial decomposition isgo, iko, io{k), k = 1, ..., 
iko, of the integer K^ - K2 by subtracting the symbolic 
polynomial representation of K2 from that of K^ (as 
polynomials) to produce a symbolic polynomial repre- 
sentation ofK^ - K2 from which the polynomial decom- 
position ofK^ - K2 can be obtained. Here nkmax is the 
dimension of the arrays /a, ib, io in the calling routine 
and mhalf equals 2^^ It is noted that by setting isgb 
equal to -isgb the polynomial decomposition ofK^ + K2 
can also be obtained with this subroutine. 

subroutine muldif(ia, ib, io, isga, isgb, isgo, ika, ikb, iko, 
* nkmax,mhalf) 

integer ia(*), ib(*), io(*) 

integer isga, isgb, isgo, ika, ikb, iko, nkmax, mhalf 
integer i, ikol, irel 
if(isgb.eq.O)then 



100 



200 



if(isga.eq.O)then 

isgo=0 

iko = 2 

io(l) = 

io(2) = 

return 
endif 

isgo = isga 
iko = ika 
do 100i=l,iko 

io(i) = ia(i) 
continue 
elseif(isga.eq.O)then 
isgo =-isgb 
iko = ikb 
do 200 i=l,iko 

io(i) = ib(i) 
continue 



else 

iko = ika 
if(ikb.lt.ika) then 

do300i=ikb+l,ika 
ib(i) = 
300 continue 

elseif(ika.lt.ikb) then 
iko = ikb 

do400i=ika+l,ikb 
ia(i) = 
400 continue 

endif 

if(isga*isgb.gt.O)then 
irel = 

do 500 i = iko, 1,-1 
if(ia(i).gt.ib(i))then 
irel = 1 
go to 600 
elseif(ia(i).h.ib(i))then 
irel = -l 
go to 600 
endif 
500 continue 

600 continue 

if(irel.eq.O)then 
isgo = 
do 700 i=l,iko 
io(i) = 
700 continue 

else 

isgo=isga*irel 

io(l) = (ia(l)-ib(l))*irel 

do 800 i=2,iko 

if(io(i-l).lt.O)then 
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800 



io(i)=-l 

io(i-l) = io(i-l) + mlialf 
else 

io(i) = 
endif 

io(i) = io(i) + (ia(i)-ib(i))*irel 
continue 
if(io(iko).lt.O) stop 210 
endif 
else 

isgo=isga 

io(l) = ia(l)+ib(l) 

do 900 i=2,iko 

if(io(i-l).ge.mhalf) then 
io(i)=l 

io(i-l) = io(i-l) - mhalf 
else 

io(i) = 
endif 

io(i) = io(i) + ia(i)+ib(i) 
900 continue 

if(io(iko).ge.nihalf) then 
iko = iko+1 

if(iko.gt.nkmax) stop 220 
io(iko) = 1 

io(iko-l) = io(iko-l) - mhalf 
endif 
endif 
endif 

if(iko .eq. 2) go to 990 
ikol =iko 
do950i = ikol,3,-l 

if(io(i) .ne. 0) go to 990 
iko = iko - 1 
950 continue 
990 continue 
retum 
end 



6. Application: Locating a Point Relative 
to a Plane 

Given an integer n,n>4, let 5* be a set ofn points in 
3 -dimensional space. Given an integer /, < / < 9, let x,, 
y^, z^, / = 1, ..., /2, be the input decimal coordinates of the 
points in S, and assume that their double precision 
floating point representations have been rounded off at 
the /th digit of their decimal fractions through applica- 
tions, one per coordinate, of the two-integer decompo- 
sition process. Accordingly, let x(i), ix(i), ix2(i), y(i), 
iy(i), iy2(i), z(i), iz(i), iz2(i), i = 1, ..., n, be the arrays 



produced by the three applications of the two-integer 
decomposition process that contain the rounded off x-, 
y-, z-coordinates and their two-integer decompositions. 

Given points pi, p^^ p^, in S that are vertices of a non- 
degenerate triangle, a fundamental problem in compu- 
tational geometry is that of finding the location of a 
point />4 in S relative to the plane //that contains the tri- 
angle. Let ¥t be the open half-space defined by H for 
whichjt?!, />2, /?3 appear in a counterclockwise direction 
around the boundary of the triangle when looking at the 
triangle from ¥t. Let //" be the other half-space defined 
by //. Determining in which of //, /T, //", the point p^ 
is located may not on occasion be satisfactorily done 
using floating point arithmetic. Accordingly, the fol- 
lowing is a Fortran subroutine called crsinn for doing 
this using polynomial decompositions. On output the 
sign isgo (-1, 0, 1) of some polynomial decomposition 
determines the location oip^ (//", //, /T). 

This routine actually does more. It produces polyno- 
mial decompositions isgox^ ikox, iox(k), k= 1, ..., ikox, 
^sgoy, ikoy, ioy(k), k = 1, ..., ikoy, isgoz, ikoz, ioz(k), 
k= i, ..., ikoz, of integers that are the coordinates of a 
vector V pointing into H^ and perpendicular to //. It also 
produces the polynomial decomposition isgo, iko, io{k), 
k= 1, ..., iko, of an integer whose sign isgo determines 
the location of P4 and whose value when divided by 
both 10^ and the length of v is the perpendicular dis- 
tance from/>4 to //. Here mhalf equals l}^, mfull equals 
2^^, and ifir, isec, ithi, ifou are locations in the arrays ix, 
ixl, etc. corresponding to the points py, p^, Pz, Pa, 
respectively. In addition, isclp{k), k= 1, 2, is an array 
such that the polynomial decomposition of 10^ is 
isclp{\), isclpil) (the sign of 10^ and the dimension of 
array isclp are already known to be 1 and 2, respective- 

ly). 

subroutine crsinn(ix, iy, iz, ix2, iy2, iz2, ifir, isec, ithi, 

* ifou, mhalf, mfull, isclp, io, isgo, iko, iox, 

* isgox, ikox, ioy, isgoy, ikoy, ioz, isgoz, ikoz) 
integer ix(*), iy(*), iz(*), ix2(*), iy2(*), iz2(*) 

integer io(*), iox(*),ioy(*), ioz(*) 
integer ifir, isec, ithi, ifou 
integer isclp(*), mhalf, mfull, nkmax 
parameter (nkmax = 30) 
integer iu(nkmax), iv(nkmax), iw(nkmax) 
integer ixt(nkmax), iyt(nkmax), izt(nkmax) 
integer ix3(nkmax), iy3(nkmax), iz3(nkmax) 
integer ix4(nkmax), iy4(nkmax), iz4(nkmax) 
integer ixf(nkmax), iyf(nkmax), izf(nkmax) 
integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew 
integer ixthw, iythw, izthw, ixfow, iyfow, izfow 
integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2 
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integer ixth2, iyth2, izth2, ixfo2, iyfo2, izfo2 

integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf 

integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2 

integer isgx3, isgy3, isgz3, ikx3, iky3, ikz3 

integer isgx4, isgy4, isgz4, ikx4, iky4, ikz4 

integer isgo, iko, isgox, ikox, isgoy, ikoy, isgoz, ikoz 

integer isgu, isgv, isgw, iku, ikv, ikw 

ixfiw = ix(ifir) 

iyfiw = iy(ifir) 

izfiw = iz(ifir) 

ixsew = ix(isec) 

iysew = iy(isec) 

izsew = iz(isec) 

ixthw = ix(ithi) 

iythw = iy(ithi) 

izthw = iz(ithi) 

ixfow = ix(ifou) 

iyfow = iy(ifou) 

izfow = iz(ifou) 

ixfi2 = ix2(ifir) 

iyfi2 = iy2(ifir) 

izfi2 = iz2(ifir) 

ixse2 = ix2(isec) 

iyse2 = iy2(isec) 

izse2 = iz2(isec) 

ixth2 = ix2(ithi) 

iyth2 = iy2(ithi) 

izth2 = iz2(ithi) 

ixfo2 = ix2(ifou) 

iyfo2 = iy2(ifou) 

izfo2 = iz2(ifoxi) 

call decmp2(ixf, isgxf, ikxf, ixfiw, ixfi2, mhalf, mfull, isclp) 

call decmp2(iyf, isgyf, ikyf, iyfiw, iyfi2, mhalf, mfull, isclp) 

call decmp2(izf, isgzf, ikzf, izfiw, izfi2, mhalf, mfull, isclp) 

call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp) 

call muldif(io, ixf, ixt, isgo, isgxf, isgx2, iko, ikxf, ikx2, 

nkmax,mhalf) 
call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp) 
call muldif(io, iyf, iyt, isgo, isgyf, isgy2, iko, ikyf, iky2, 

nkmax,mhalf) 
call decmp2(io, isgo, iko, izsew, izse2, mhalf, mfull, isclp) 
call muldif(io, izf, izt, isgo, isgzf, isgz2, iko, ikzf, ikz2, 

nkmax,mhalf) 
call decmp2(io, isgo, iko, ixthw, ixtli2, mhalf, mfull, isclp) 
call muldif(io, ixf, ix3, isgo, isgxf, isgx3, iko, ikxf, ikx3, 

nkmax,mhalf) 
call decmp2(io, isgo, iko, iythw, iyth2, mhalf, mfull, isclp) 
call muldif(io, iyf, iy3, isgo, isgyf, isgy3, iko, ikyf, iky3, 

nkmax,mhalf) 
call decmp2(io, isgo, iko, izthw, izth2, mhalf, mfull, isclp) 
call muldif(io, izf, iz3, isgo, isgzf, isgz3, iko, ikzf, ikz3, 

nkmax,mhalf) 



call decmp2(io, isgo, iko, ixfow, ixfo2, mhalf, mfull, isclp) 
call muldif(io, ixf, ix4, isgo, isgxf, isgx4, iko, ikxf, ikx4, 

* nkmax,mhalf) 

call decmp2(io, isgo, iko, iyfow, iyfo2, mhalf, mfull, isclp) 
call muldif(io, iyf, iy4, isgo, isgyf, isgy4, iko, ikyf, iky4, 

* nkmax,mhalf) 

call decmp2(io, isgo, iko, izfow, izfo2, mhalf, mfull, isclp) 
call muldif(io, izf, iz4, isgo, isgzf, isgz4, iko, ikzf, ikz4, 

* nkmax,mhalf) 

call mulmul(iyt, iz3, iv, isgy2, isgz3, isgv, iky2, ikz3, ikv, 

* nkmax,mhalf) 

call mulmul(izt, iy3, iu, isgz2, isgy3, isgu, ikz2, iky3, iku, 

* nkmax,mhalf) 

call muldif(iv, iu, iox, isgv, isgu, isgox, ikv, iku, ikox, 

* nkmax,mhalf) 

call mulmul(iox, ix4, io, isgox, isgx4, isgo, ikox, ikx4, iko, 

* nkmax,mhalf) 

call mulmul(izt, ix3, iv, isgz2, isgx3, isgv, ikz2, ikx3, ikv, 

* nkmax,mhalf) 

call mulmul(ixt, iz3, iu, isgx2, isgz3, isgu, ikx2, ikz3, iku, 

* nkmax,mhalf) 

call muldif(iv, iu, ioy, isgv, isgu, isgoy, ikv, iku, ikoy, 

* nkmax,mhalf) 

call mulmul(ioy, iy4, iu, isgoy, isgy4, isgu, ikoy, iky4, iku, 

* nkmax,mhalf) 
isgu =-isgu 

call muldif(io, iu, iw, isgo, isgu, isgw, iko, iku, ikw, 

* nkmax,mhalf) 

call mulmul(ixt, iy3, iv, isgx2, isgy3, isgv, ikx2, iky3, ikv, 

* nkmax,mhalf) 

call mulmul(iyt, ix3, iu, isgy2, isgx3, isgu, iky2, ikx3, iku, 

* nkmax,mhalf) 

call muldif(iv, iu, ioz, isgv, isgu, isgoz, ikv, iku, ikoz, 

* nkmax,mhalf) 

call mulmul(ioz, iz4, iu, isgoz, isgz4, isgu, ikoz, ikz4, iku, 

* nkmax,mhalf) 
isgu =-isgu 

call muldif(iw, iu, io, isgw, isgu, isgo, ikw, iku, iko, 

* nkmax,mhalf) 
retum 

end 

Sometimes besides knowing the location of the point 
jt?4 relative to the plane H it may be desirable to know 
the perpendicular distance from p^ to H. The following 
is Fortran code for this purpose. It uses the polynomial 
decompositions that are part of the output of subroutine 
crsinn. Variables here are either integer or double pre- 
cision following convention. Here r215 equals 
{2MQy\ dscle equals {\QMQ)\ and dist is the result- 
ing signed perpendicular distance. In addition, it is 
assumed that subroutine doubnm (presented below) 
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exists for transforming the polynomial decomposition 
of an integer into the double precision floating point 
value of the integer. 



twelve points are read into double precision arrays x(i), 
y(i), z(i), / = 1, ..., 12, so that x(i), y(i), z(i) contain the 
x-,y-, z-coordinates, respectively, of point z. 



call crsinn(ix, iy, iz, ix2, iy2, iz2, ifir, isec, ithi, ifou, 

* mhalf, mfUll, isclp, io, isgo, iko, iox, isgox, 

* ikox, ioy, isgoy, ikoy, ioz, isgoz, ikoz) 
call doubnm(io, isgo, iko, r215, dnum) 

call doubnm(iox, isgox, ikox, r215, xnum) 

call doubnm(ioy, isgoy, ikoy, r215, ynum) 

call doubnm(ioz, isgoz, ikoz, r215, znum) 

dnux = dmaxl(dabs(xnum),dabs(ynum),dabs(znum)) 

xnum = xnum/dnux 

ynum = ynum/dnux 

znum = znum/dnux 

dnom = dsqrt(xnum**2+ynum**2+znum**2) 

dist = ((dnum/dnux)/dnom)/dscle 

The following is subroutine doubnm that was called 
above. 

subroutine doubnm(io, isgo, iko, r215, dnum) 

integer io(*) 

double precision dnum, r215, rpwr 

integer isgo, iko, i 

if(isgo.eq.O) then 

dnum = O.OdO 

go to 900 
else 

if(iko.lt. 2)stop310 

if(iko .gt. 68) stop 320 

rpwr= l.OdO 

dnum = dble(io(l)) 

do 100 i = 2, iko 
rpwr = rpwr*r215 
dnum = dnum + dble(io(i))*rpwr 
100 continue 
endif 

if(isgo.lt.O) dnum = -dnum 
900 continue 
return 
end 



7. Numerical Examples 

Twelve lines follow, each line containing three num- 
bers. Each line corresponds to a point in 3-dimensional 
space, and the three numbers in the line correspond to 
the X-, y-, z-coordinates of the point, in that order. 
Given z, 1 < z < 12, point / is the point corresponding to 
the /th line. It is assumed that the coordinates of the 



-13.729277089 

38.000000000 

41.736468803 

85.557213025 

33.675274550 

1.724283838 

15.161728368 

0.082570927 

47.541325082 

-33.285508962 

-2.277195916 

70.061142979 



14.530621914 

7.049967880 

68.831641719 

^9.840807038 

-77.937397763 

-53.594476834 

3.186043237 

-30.956721161 

-77.446759923 

-14.545102894 

-58.886394970 

9.068097315 



97.981467003 
-92.123710427 
-59.331882431 
-13.994897166 

52.741164465 
-84.424190762 

98.792566086 
-95.085758310 
^1.735139045 

93.175307798 

80.791131020 
-70.800333278 



Given / equal to 8, through the two-integer decompo- 
sition process, the numbers above are rounded off at the 
/th = 8th digit of their decimal fractions and saved in 
x(i), y(i), z(i), i = 1, ..., 12, so that then they appear as 
follows. 



-13.729277090 
38.000000000 
41.736468800 
85.557213020 
33.675274550 
1.724283840 
15.161728370 
0.082570930 
47.541325080 

-33.285508960 
-2.277195920 
70.061142980 



14.530621910 

7.049967880 

68.831641720 

^9.840807040 

-77.937397760 

-53.594476830 

3.186043240 

-30.956721160 

-77.446759920 

-14.545102890 

-58.886394970 

9.068097310 



97.981467000 
-92.123710430 
-59.331882430 
-13.994897170 

52.741164460 
-84.424190760 

98.792566090 
-95.085758310 
^1.735139040 

93.175307800 

80.791131020 
-70.800333280 



Each rounded off coordinate is also decomposed into 
two integers. Twelve lines follow, each line containing 
six integers. For each /, i = 1, ..., 12, the first two inte- 
gers in the /th line are the two integers into which the x- 
coordinate of point / is decomposed. Similarly the next 
two integers correspond to the j^-coordinate, and the 
last two to the z-coordinate. The two-integer decompo- 
sitions of the twelve points are then saved into /x(/), 
ix2(i), iy(i), iy2(i), iz(i), zz2(/), i = 1, ..., 12, in the obvi- 
ous manner. It is noted that when mfull = 2^^ = 
1073741824 appears as the second integer correspon- 
ding to a coordinate it is to be interpreted as a zero. 

-13 -72927709 14 53062191 97 98146700 
38 1073741824 704996788 -92 -12371043 
41 73646880 68 83164172 -59 -33188243 
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85 


55721302 


^9 


-84080704 


-13 


-99489717 


33 


67527455 


-77 


-93739776 


52 


74116446 


172428384 





-53 


-59447683 


-84 ^2419076 


15 


16172837 


318604324 





98 


79256609 


8257093 





-30 


-95672116 


-95 


-8575831 


47 


54132508 


-11 ^4675992 


-A\ 


-73513904 


-33 


-28550896 


-14 


-54510289 


93 


17530780 


227719592 





-58 


-88639497 


80 


79113102 


70 


6114298 


906809731 





-70 


-80033328 



Given / as above equal to 8 and setting isclu to 10^ = 
10^, by calling subroutine decomp the polynomial 
decomposition of isclu is found to be isgcl, isclp{\), 
isclp(2), with isgcl equal to 1, isclp{\) equal to 24832 
and isclp{2) equal to 3051 (10^ = 24832 + 3051 ■ mhalf 
= 24832 + 3051 ■ 2^^). 

As an example of how exact computations are car- 
ried out that reference coordinates of the input points 
given above, the computation that is the product of the 
x-coordinate of point 2 and the j-coordinate of point 4 
minus the product of the j-coordinate of point 2 and the 
x-coordinate of point 4 is described. Using rounded-off 
numbers the result of the computation should equal 
x(2)-j(4)-j(2)-x(4),i.e., 

(38.000000000) ■ H9.840807040) - (7.049967880) ■ (85.557213020). 

On the other hand, if exact computations are required 
then each of the four numbers involved must first be 
converted into an integer which is the number times 1 0^ 
rounded off at the decimal point. Since the resulting 
integer may be too big to be saved into a 32 bit integer 
word, its polynomial decomposition in the form of two 
or more 32 bit integer words is obtained instead. Thus, 
for example, the polynomial decomposition of the inte- 
ger which is the x-coordinate of point 4 times 10^ 
rounded off at the decimal point can be obtained by 
calling subroutine decmp2 using the two-integer 
decomposition of the coordinate, i.e. the integers 85 
and 55721302, and the polynomial decomposition of 
10^ as obtained above. The resulting polynomial 
decomposition is then found to be isgoxA, ikoxA, 
iox4(k), k= 1, ..., ikox4, with isgox4 equal to 1, ikox4 
equal to 3, iox4(k), ^ = 1, 2, 3, equal to 29270, 31723, 
and 7 (the integer which is the x-coordinate of point 4 
times 10^ rounded off at the decimal point equals 
29270 + 31723 ■ (2^^) + 7 ■ (2^^). In the same manner 
the polynomial decompositions associated with the x-, 
^-coordinates of point 2, and the j-coordinate of point 
4 are found to be, respectively, isgoxl, ikoxl, iox2(k), 
k= 1, ..., /^ox2, isgoyl, ikoyl, ioy2{k), k= I, ..., ikoy2, 
isgoy4, ikoy4, ioy4(k), k = 1, ..., ikoy4, with isgox2 



equal to 1, /^ox2 equal to 3, iox2(k), ^ = 1, 2, 3, equal 
to 26112, 17662, 3, isgoy2 equal to 1, ikoy2 equal to 2, 
ioy2(k),k= 1,2, equal to 26036, 21514, isgoy4 equal to 
-1, ikoy4 equal to 3, ioy4(k), k= 1,2, 3, equal to 2368, 
21030, 4. Finally, using these polynomial decomposi- 
tions as input the desired result is obtained by calling 
subroutines mulmul and muldif as follows. Here nkmax 
is the dimension of all of the arrays (input and output). 

callmulmul(iox2,ioy4,iu,isgox2,isgoy4,isgu,ikox24koy4,iku,nkmax,mhalf) 
callmulmul(ioy2,iox4,iv,isgoy2,isgox4,isgv,ikoy2,ikox4,ikv,nkmax,mhalf) 
callmuldif(iu,iv,io,isgu,isgv,isgo,iku,ikv,iko,nkmax,mhalf) 

The polynomial decomposition of an integer that 
approximates the desired result times (10^)^ rounded off 
at the decimal point is then found to be isgo, iko, io(k), 
k=l, ..., iko, with isgo equal to -1 , iko equal to 5, io{k), 
k=\, ..., 5, equal to 21112, 15183,31880,21597,21. 
The desired result is then approximately equal to this 
integer, i.e. 21112 + 15183 ■ (2^0 + 31880 ■ (2^0' + 
21597 ■ (2^^)^ + 21 ■ (2'Y, divided by (10^)'. By calling 
subroutine doubnm the integer can be approximated by 
a double precision number which when divided by 
(10^)^ is approximately -2497. 12627 12133, an approx- 
imation of the desired result. 

Finally, an example is given for the purpose of 
describing the process of locating a point relative to a 
plane that contains three other points. Here the point 
whose location is desired is point 12 as given above, 
and the three other points are point 1 , point 2, point 8 
also as given above. With / as the triangle with vertices 
point 1, point 2, point 8, and H as the plane that con- 
tains /, H^ is taken to be the open half-space defined by 
H for which point 1 , point 2, point 8 appear in a coun- 
terclockwise direction around the boundary of t when 
looking at t from H^. /T is taken to be the other half- 
space defined by H. With ifir set to 1 , isec set to 2, ithi 
set to 8, ifou set to 12, which of 77, H^, H~ contains point 
12 can be determined by calling subroutine crsinn as 
follows. 

call crsimi(ix, iy, iz, ix2, iy2, iz2, ifir, isec, ithi, ifou, 

* mhalf, mfUll, isclp, io, isgo, iko, iox, isgox, 

* ikox, ioy, isgoy, ikoy, ioz, isgoz, ikoz) 

On output isgo equals -1 so that point 12 must be in H~. 
In addition iko equals 7, io{k\ k = \, ..., 7, equals 
21844, 27853, 3870, 5372, 13630, 9887, 213, isgox 
equals -1, ikox equals 5, iox{k), k = \, ..., 5, equals 
11868, 15341, 2677, 15631, 62, isgoy equals 1, ikoy 
equals 5, ioy{k), k=\, ..., 5, equals 11577, 8756, 364, 
27887, 63, isgoz equals -1, ikoz equals 5, ioz{k), k=\, 
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..., 5, equals 5921, 23919, 26934, 16812, 19. By calling 
subroutine doubnm the integer whose polynomial 
decomposition is isgo, iko, io(k), k= 1, ..., iko, can be 
approximated by a double precision number dnum. 
Similarly, the three integers, say zx, iy, iz, whose poly- 
nomial decompositions are isgox, ikox, iox(k), k= 1, ..., 
ikox, isgoy, ikoy, ioy(k), k = 1, ..., ikoy, isgoz, ikoz, 
ioz{k), k= 1, ..., ikoz, can be approximated, respective- 
ly, by double precision numbers xnum, ynum, znum. 
The vector (/x, iy, iz), which is perpendicular to the 
plane H and points into H^, can then be approximated 
by the vector {xnum, ynum, znum). Finally by dividing 
dnum by both the length of the vector (xnum, ynum, 
znum) and 10^ the number -25.047402554921 is 
approximately obtained, an approximation of the 
signed perpendicular distance from point 1 2 to plane //, 
the negative sign indicating that point 1 2 is in /T. 
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8. Summary 

A scheme has been presented and software has been 
documented for transforming into a series of integers 
input decimal numbers that have been read into a com- 
puter as double precision floating point numbers and 
for carrying out multiplication, addition and subtraction 
operations based on these numbers using their integer 
representations. The total number of significant digits 
of each input number must not be greater than 14, and 
the number of digits to the left of the decimal point 
must not exceed 9. Through a preprocessing step the 
double precision floating point representation of each 
input decimal number is rounded off at a prespecified 
location of its decimal fraction, a location no more than 
9 digits to the right of the decimal point, and the round- 
ed off number is decomposed into two integers. All 
operations that follow involving the input number are 
then carried out in terms of the rounded off double pre- 
cision floating point number and when this is not satis- 
factory in terms of the two integers. This scheme has 
been successfully incorporated into several computa- 
tional geometry related programs such as REGTET [1]. 
Other programs that incorporate this scheme can be 
found at http://math.nist.govrJBemal/JBemal_Sft.html. 
Besides being used for locating a point relative to a 
plane, in these programs the scheme has also been used 
for locating a point relative to a sphere, for computing 
the intersection of a line and a plane, for computing the 
center of a sphere, etc. 
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